% ***************************************************************************************
% m_identificationtest.m 
% ***************************************************************************************
clear all; 
close all;
format shortG;

% ==================== SET FUNDAMENTAL PARAMETERS ========================
load parameters.mat;                           
% ========================= LOAD DATA  =========================================  
load calib_setup.mat;                    
N = size(Rdata_all,1); 
NT = size(Rdata_all,2)-1;   
Rdensdata = Rdata_all(:,2:end)./repmat(Area,1,NT);  
Ldensdata = Ldata_all(:,2:end)./repmat(Area,1,NT); 

% ================= LOAD COMMUTING COST ===================================
load commuting_cost.mat;  
% ============= LOAD INVERSION RESULTS FROM MAIN CALIBRATION ==============
load calib_auxiliary.mat;    
Xmat = [zeros(N,1), X55, X60, X65, X70, X75]; 
Ymat = [zeros(N,1), Y55, Y60, Y65, Y70, Y75]; 

% =========== (i) DECOMPOSE FUNDAMENTALS FOR ESTIMATED VALUES ===============
% a0: log of fundamental productivity matrix (N x NT-1) from 1955,60,65,70,75
% a0_err: a0 minus time trend of a0 (log scale) N x NT-1
% b0: log of fundamental amenity matrix (N x NT-1) from 1955,60,65,70,75
% b0_err: b0 minus time trend of b0 (log scale) N x NT-1
[a0,a0_err] = ffund2(alphaest,Ymat,Ldensdata,rho,sigma,thetavec,N,NT);     
[b0,b0_err] = ffund2(betaest,Xmat,Rdensdata,rho,sigma,thetavec,N,NT);
d_a0_err = a0_err(:,2:NT-1)-a0_err(:,1:NT-2);  % N x NT-2 
d_b0_err = b0_err(:,2:NT-1)-b0_err(:,1:NT-2);
d_a0 = sum(d_a0_err,2)./(NT-2); 
d_b0 = sum(d_b0_err,2)./(NT-2);
a0 = sum(a0,2)./(NT-1); 
b0 = sum(b0,2)./(NT-1); 
% =========== (ii) DECOMPOSE FUNDAMENTALS FOR WEAK VALUES  ================
[a1,a1_err] = ffund2(0.01,Ymat,Ldensdata,rho,sigma,thetavec,N,NT);     
[b1,b1_err] = ffund2(0.01,Xmat,Rdensdata,rho,sigma,thetavec,N,NT);
d_a1_err = a1_err(:,2:NT-1)-a1_err(:,1:NT-2);   
d_b1_err = b1_err(:,2:NT-1)-b1_err(:,1:NT-2);
d_a1 = sum(d_a1_err,2)./(NT-2); 
d_b1 = sum(d_b1_err,2)./(NT-2);
a1 = sum(a1,2)./(NT-1); 
b1 = sum(b1,2)./(NT-1); 

% =========== (iii) DECOMPOSE FUNDAMENTALS FOR STRONG VALUES  ================
[a2,a2_err] = ffund2(0.3,Ymat,Ldensdata,rho,sigma,thetavec,N,NT);     
[b2,b2_err] = ffund2(0.3,Xmat,Rdensdata,rho,sigma,thetavec,N,NT);
d_a2_err = a2_err(:,2:NT-1)-a2_err(:,1:NT-2);   
d_b2_err = b2_err(:,2:NT-1)-b2_err(:,1:NT-2);
d_a2 = sum(d_a2_err,2)./(NT-2); 
d_b2 = sum(d_b2_err,2)./(NT-2);
a2 = sum(a2,2)./(NT-1); 
b2 = sum(b2,2)./(NT-1); 

% =========== FIGURES =====================================================
d_b_mat = [d_b0, d_b1, d_b2]; 
d_a_mat = [d_a0, d_a1, d_a2];

resultmat=array2table([cbddist, Area, d_b_mat, d_a_mat],'VariableNames', {'cbddist', 'area','db_base', 'db_weak','db_strong','da_base','da_weak','da_strong'}); 
writetable(resultmat,'../../output/quant/output_calib_identidicationtest_change.csv');

b_mat = [b0, b1, b2]; 
a_mat = [a0, a1, a2];

resultmat=array2table([cbddist, Area, b_mat, a_mat],'VariableNames', {'cbddist', 'area','b_base', 'b_weak','b_strong','a_base','a_weak','a_strong'}); 
writetable(resultmat,'../../output/quant/output_calib_identidicationtest_level.csv');


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%                       FUNCTIONS
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% --------  FUNCTION FOR COMPUTING THE EXOGENOUS FUNDAMENTALS --------
function[c, c_err] = ffund2(param,Zmat,densdata,rho,sigma,thetavec,N,NT)  

Z0=Zmat(:,NT); cmat=zeros(N,NT-1);
for t=1:NT-1
    denspost=densdata(:,NT-t+1); 
    invspill = denspost.^(-param);
    cc = invspill.*Z0;
    cmat(:,NT-t)= cc; 
    theta=thetavec(NT-t); 
    Z0=Zmat(:,NT-t).*(Z0.^(-rho*(1-theta)));     
end

c = log(cmat);  
c_te = sum(c,1)./N;
c_err = c-repmat(c_te,N,1); 
end 

